Atoms of multistationarity in chemical reaction networks 
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Abstract 

Chemical reaction systems are dynamical systems that arise in chemical engineering and 
systems biology. In this work, we consider the question of whether the minimal (in a pre- 
cise sense) multistationary chemical reaction networks, which we propose to call 'atoms of 
multistationarity,' characterize the entire set of multistationary networks. Our main result 
states that the answer to this question is 'yes' in the context of fully open continuous-flow 
stirred-tank reactors (CFSTRs), which are networks in which all chemical species take part 
in the inflow and outflow. In order to prove this result, we show that if a subnetwork 
admits multiple steady states, then these steady states can be lifted to a larger network, 
provided that the two networks share the same stoichiometric subspace. We also prove an 
analogous result when a smaller network is obtained from a larger network by 'removing 
species.' Our results provide the mathematical foundation for a technique used by Siegal- 
Gaskins et al. of establishing bistability by way of 'network ancestry' Additionally, our 
work provides sufficient conditions for establishing multistationarity by way of atoms and 
moreover reduces the problem of classifying multistationary CFSTRs to that of cataloging 
atoms of multistationarity. As an application, we enumerate and classify all 386 bimolecular 
and reversible two-reaction networks. Of these, exactly 35 admit multiple positive steady 
states. Moreover, each admits a unique minimal multistationary subnetwork, and these 
subnetworks form a poset (with respect to the relation of 'removing species') which has 11 
minimal elements (the atoms of multistationarity). 

Keywords: chemical reaction networks, mass-action kinetics, multiple steady states, 
Jacobian Criterion, injectivity 



1 Introduction 

This work concerns an important class of dy- 
namical systems arising in chemical engineer- 
ing and systems biology, namely, chemical re- 
action systems. As bistable chemical systems 
are thought to be the underpinnings of bio- 
chemical switches, a key question is to de- 
termine which systems admit multiple steady 
states. In this work, we consider the question 




Figure 1: We propose to call 2B — s- A — >• A + B an 
'atom of multistationarity'; see Section [5] 
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of whether the minimal (in a precise sense) networks, which we propose to call 'atoms of multi- 
stationarity,' characterize the entire set of multistationary networks. We prove that such atoms 
do characterize multistationarity for the case of fully open continuous-flow stirred-tank reac- 
tors (CFSTRs), which are networks in which all chemical species take part in the inflow and 
outflow (see Definition 2.4). For instance, the five networks depicted in Figure [T] are multista- 
tionary in the CFSTR setting, but only one is minimal with respect to 'removing species' (see 
Theorem 7.1). Following other analyses of small networks \1'6\ [2"T] [2"2"1 |2"3"1 [2"1"1 [2"5] , our main 
application is to those CFSTRs containing two non-flow reactions. 

Chemical reaction systems are nonlinear and parametrized by unknown reaction rate con- 
stants. Thus, determining whether a chemical reaction network admits multiple steady states 
is difficult: for instance, in the mass-action kinetics setting, it requires determining existence 
of multiple positive solutions to a system of polynomials with unknown coefficients. However, 
various criteria have been developed that often can answer this question. For instance, the 
Deficiency, Advanced Deficiency, and Higher Deficiency Theories developed by Ellison, Fein- 
berg, Horn, Jackson, and Ji in many cases can affirm that a network admits multiple steady 
states or can rule out the possibility [9j HU El H5] . Similarly, the Jacobian Criterion and the 
more general injectivity test developed by Craciun and Feinberg can preclude multiple steady 
states 01 El El ESI EZj- These results have been implemented in the CRN Toolbox, freely 
available computer software developed by Feinberg and improved by Ellison and Ji |10| . Re- 
lated software programs include BioNetX |19l [2UJ and Chemical Reaction Network Software for 
Mathematica |18j . 

For systems for which the above software approaches are inconclusive, Conradi et al. advo- 
cate an approach which first determines whether certain subnetworks admit multiple positive 
steady states, and if so, tests whether these instances can be lifted to the original network [2]. 
Here, we too examine the topic of lifting multistationarity from a subnetwork to an overall 
network: Theorem 3A_ states that this can be accomplished as long as the steady states of 
interest are nondegenerate and the two networks share the same stoichiometric subspace. This 
result and its proof extend Theorem 2 in work of Craciun and Feinberg [3] . An important con- 
sequence of our theorem is that it provides the mathematical foundation for the technique of 
Siegal-Gaskins et al. which establishes bistability by way of 'network ancestry' (see Remark 3.3 ); 
their method was applied to a large class of simple gene regulatory networks |23j . 

A fully open continuous-flow stirred-tank reactor (CFSTR) is a network in which all chem- 
ical species enter the system at constant rates and are removed at rates proportional to their 
concentrations (see Definition 2.4). In the setting of these systems, we extend our results be- 
yond subnetworks to 'embedded networks' which are obtained by removing species as well as 
reactions from a network (Definition 2.2). Corollary 4.6 states that if an embedded CFSTR of 
a CFSTR is multistationary then so is the CFSTR itself. Therefore, the set of multistation- 
ary CFSTRs is characterized by its minimal elements (with respect to the embedded network 
relation), and we pose the challenge of characterizing these atoms. In this work, we focus on 
cataloging the smallest atoms. 

Recent work of the first author presented a simple characterization of the one-reaction fully 
open CFSTRs that admit multiple steady states in the mass-action kinetics setting (Theo- 
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2.11) [16]. Here we consider the bimolecular two-reaction CFSTRs; a network is 'bimolec- 



rem 

ular' if all of its chemical complexes contain at most two molecules. We enumerate all 386 
reversible such networks. Of these, exactly 35 admit multiple positive steady states. Moreover, 
each admits a unique minimal multistationary sub-CFSTR, and these subnetworks form a poset 
with respect to the embedded network relation that has 11 minimal elements (Theorem 7.1). 
These 11 networks are precisely the CFSTR atoms of multistationarity in the bimolecular two- 
reaction setting. Note that a similar enumeration of small bimolecular networks was undertaken 
by Deckard, Bergmann, and Sauro jS], from which Pantea and Craciun sampled networks to 
compute the fraction of such networks that pass the Jacobian Criterion [20J. 

This article is organized as follows. Section [2] introduces chemical reaction systems. Our 
main result for lifting multiple steady states from subnetworks, Theorem 3.1, appears in Sec- 
tion [3} In Section [4j this result is extended in the case of fully open CFSTRs: Theorem 4.2 



implies that steady states from embedded CFSTRs can be lifted as well. Section [5] introduces 
'atoms of multistationarity' Section [6] describes our approach to enumerating bimolecular 
two-reaction networks (Algorithm |6.4| ), and Section [7] determines which such networks are mul- 
tistationary in the mass-action kinetics setting (Theorem 7.1) and displays the resulting atoms 
of multistationarity (Figure f3j. 



2 Chemical reaction network theory 

In this section we review the standard notation and recall the classification of one-reaction 
CFSTRs. 



2.1 Chemical reaction networks 

We begin with an example of a chemical reaction: 

2X X +X 2 -> X 3 . 

Each Xi is called a chemical species, and 2Xi + X2 and X3 are called chemical complexes. 
Assigning the reactant complex 2X\ + X2 to the vector y = (2, 1,0) and the product complex 
X3 to the vector y' = (0, 0, 1), we can write the reaction as y — > y' . In general we let s denote 
the total number of species Xi, and we consider a set of r reactions, each denoted by 

Vk -»• y' k , 

for k £ {1, 2, . . . , r}, and yk,y'f. £ ^>o> with yt 7^ y' k - We index the entries of a complex vector 
Vk by writing y^ = (yki,Vk2, ■ ■ ■ iVks) £ ^>o> anc ^ we wm can Vki the stoichiometric coefficient 
of species i in complex y k . For ease of notation, when there is no need for enumeration we 
typically will drop the subscript k from the notation for the complexes and reactions. 

Definition 2.1. Let S = {Xi}, C = {y}, and 1Z = {y — > y'} denote finite sets of species, 
complexes, and reactions, respectively. The triple {S, C, 1Z} is called a chemical reaction network 
if it satisfies the following: 
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1. for each complex y G C, there exists a reaction in TZ for which y is the reactant complex 
or y is the product complex, and 

2. for each species Xi £ S, there exists a complex y £ C that contains Xi. 

A network decouples if there exist nonempty subsets TZ' C TZ and TZ" C TZ such that 
TZ = TZ' U TZ" and such that the species involved in reactions in TZ' are distinct from those of 
TZ". We next define a subnetwork and the more general concept of an 'embedded' network,' 
which was introduced by the authors in |17l §4.2]. Informally, a network N is an embedded 
network of a network G if N may be obtained from G by removing reactions and 'removing 
species.' 

Definition 2.2. Let (5 = {S,C,TZ} be a chemical reaction network. 

1. Consider a subset of the species S C S, a subset of the complexes C C C, and a subset of 
the reactions R C TZ. 

• The restriction of R to S, denoted by R\s, is the set of reactions obtained by taking 
the reactions in R and removing all species not in S from the reactant and product 
complexes. If a trivial reaction (one in which the reactant and product complexes 
are the same) is obtained in this process, then it is removed. Also removed are extra 
copies of repeated reactions. 

• The restriction of C to R, denoted by C\r, is the set of (reactant and product) 
complexes of the reactions in R. 

• The restriction of S to C, denoted by S\c, is the set of species that are in the 
complexes in C. 

2. The network obtained from (3 by removing a subset of species {Xi} C S is the network 

{s\{Xi}, C\ n[sX{Xi} , TZ\ s \{ Xi }} ■ 

3. A subset of the reactions TZ' C TZ defines the subnetwork {S\c\ K/ ,C\n' ,TZ'} . 

4- Let & = {S,C,TZ} be a chemical reaction network. An embedded network of &, which 
is defined by a subset of the species, S = {X^, JQ 2 , . . . ,Xi k } C S, and a subset of the 
reactions, R = {Rj 1} Rj 2 , . . . ,Rj t } C TZ, that involve all species of S, is the network 
(S,C\r\ s , R\s) consisting of the reactions R\s- 

Remark 2.3. We note that a network is also a subnetwork and an embedded network of itself. 
In fact, any subnetwork {5|c| TC , > TZ'} is an embedded network, namely the one defined by 
the subset of species , and the subset of reactions TZ' . 

We also note for readers who are familiar with species-reaction (SR) graphs that the defi- 
nitions of 'subnetwork' and 'embedded network' can be interpreted as follows. Recall that the 
SR graph of a network consists of species vertices and reaction vertices, with edges arising from 
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reactions in the network; for details, see [5j. A subnetwork corresponds to the subgraph of the 
SR graph induced by the full set of species and the subset of reaction vertices arising from re- 
actions in the subnetwork. As for an embedded network, this arises as the subgraph induced by 
the corresponding subsets of species and reaction nodes. 

One focus of our work is on CFSTRs, which we now define. 

Definition 2.4. 1. A flow reaction contains only one molecule; such a reaction is either an 
inflow reaction — > Xi or an outflow reaction X{ — > 0. 

2. A chemical reaction network is a continuous-flow stirred-tank reactor (CFSTR) if it con- 
tains all outflow reactions X{ — )■ (for all Xi G S). A CFSTR is fully open if it addi- 
tionally contains all inflow reactions — > Xi. A sub-CFSTR is a subnetwork that is also 
a CFSTR. 

We note that Craciun and Feinberg use the term 'feed reactions' for inflow reactions and 
'true reactions' for non-flow reactions. In chemical engineering, a CFSTR refers to a well-mixed 
tank in which reactions occur. An inflow reaction represents the flow of species (at a constant 
rate) into the tank in which the non-flow reactions take place, and an outflow reaction represents 
the removal or degradation of a species (at rate proportional to its concentration). 

Example 2.5. Consider the following fully open CFSTR: 

K\ K3 K5 K7 Kg 

0t±A O^B O^C 2A^A + B^A + C. (1) 

K2 K4 KQ Kg Kio 

The following sub-CFSTR arises by removing two reactions: 

Kl K3 K 5 

+± A f± B O^C 2A^A + B^A + C. (2) 

K 2 = l K4=l K 6 = l K8 Kl ° 

Next, we obtain the following embedded network by removing species C: 

Kl K 3 

+± B 2A^ A + B <- A . (3) 

K 2 = l K 4 = l K S «X0 

2.2 Dynamics and steady states 

The concentration vector 

x(t) = (x 1 (t),x 2 {t), . . .,X\ S \{t)) 

will track the concentration Xi(t) of the i-th species at time t. A chemical reaction network 
defines a dynamical system by way of a rate function for each reaction. In other words, to each 
reaction — > y' k we assign a smooth function Rk(-) = Ry k -+y' (0 that satisfies the following 
assumption. 

Assumption 2.6. For keK, R k (-) = Ry h -+u> (•) : ^> -> K satisfies: 
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1. Ry k ^ y ' k (-) depends explicitly upon X{ only if yki 7^ 0. 

d 

2. — — Ry k ^ y ' k (x) > for those x% for which y k i 7^ ; and equality can hold only if at least 

OX{ 

one coordinate of x is zero. 

3. Ry k ^ y ' k (x) = if Xi = for some i with y k i / 0. 

A. If 1 < j/fcj < yn, then lim = where all other Xj > are held fixed in the limit. 

Xi^o R k (x) 



The final assumption simply states that if the Z-th reaction demands strictly more molecules 
of species Xi as inputs than does the k-th reaction, then the rate of the Z-th reaction decreases 
to zero faster than the fc-th reaction, as x^ — > 0. The functions R k are called the kinetics of the 
system. 

Definition 2.7. Consider a chemical reaction network {S,C,1Z = {yt —> 2/1}} and a choice of 



kinetics {R k } that satisfy Assumption 2a 



1. The following system of ODEs defines a dynamical system is called a chemical reaction 
system: 

n 

x(t) = YtRxWWk-Vk) =■ f{x(t)), (4) 
fc=i 

where the second equality is a definition. 

2. The stoichiometric subspace of the network is the span of all reaction vectors y' k — yk- We 
will denote this space by S and its dimension by a: 



S := span {y[ - y 1} 2/2 - 2/2, • • • , y'\ n \ ~ V\n\} C 



Note that Q implies that a trajectory x(t) that begins at a positive vector x(0) = c° G 
remains in the stoichiometric compatibility class, which we denote by 



OS 

So 



V := (c° + 5)nKg, (5) 

for all positive time; in other words, this set V is forward-invariant with respect to Q . 
Two points in the same stoichiometric compatibility class V are said to be stoichiometri- 
cally compatible. 

3. A concentration vector x G Mtf is a (positive) steady state of the system Q if fix) = 0. 
A steady state x is nondegenerate iflmdf(x) = S. (Here, "df(x)" is the Jacobian matrix 
of f at x: the \S\ x |<S| -matrix whose (i,j)-th entry is equal to the partial derivative 
)■ A nondegenerate steady state x is exponentially stable if each of the a := dim S 
nonzero eigenvalues of df(x) (viewed over the complex numbers) has negative real part. 
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In the case of a CFSTR, the reaction vector for the i-th inflow reaction is the i-th canonical 
basis vector of M) s \, so the stoichiometric subspace is S = M) s \. It follows that for a CFSTR, 
the unique stoichiometric compatibility class is the nonnegative orthant: V = M>q. 

An important example of kinetics is mass-action kinetics; a chemical reaction system is said 
to have mass-action kinetics if all rate functions R k take the following multiplicative form: 

R k (x) = K k xf"xf 2 ---x y { ^ =: K k x^ , (6) 
for some vector of positive reaction rate constants {k\, K2, ■ ■ ■ , i^\u\) G ^>o > W1 th the convention 



that 0° = 1. It is easily verified that each R k defined via ^ satisfies Assumption 2.6 Combining 
@ and §) gives the following system of mass-action ODEs: 

\n\ 

x(t) = 52Kkx(t)v*(j/ k -y k ) =: f(x(t)). (7) 
fc=i 

In the following example and all others in this work, we will label species by distinct letters 
such as A, B, . . . rather than X±,X2, .... 



Example 2.8. We now return to the CFSTR ([T]) in Example 2.5. The mass-action differential 
equations ([7]) for this network are the following: 



dXA 9 

= Kl — K2XA — KlX A + KgXAXB 



dt 

dXB ■> 

= K 3 - K4,X B + K 7 X A - KgXAXB - KgXAXB + K W XAXC (8) 



dt 

dxg 

dt 



K5 ~ KeX C + k 9 x a xb - n\ox A xc 



Note that a chemical reaction network gives rise to a family of mass-action kinetics systems 
parametrized by a choice of one reaction rate constant K k G M>o for each reaction, and all 
reactions not in the network can be viewed as having reaction rate constant equal to zero. We 
now generalize this concept of a parametrized family for other kinetics. 

Definition 2.9. 1. A parametrized family of kinetics K, for chemical reaction networks on 
\S\ species is an assignment to each possible reaction y k — > y' k (that involves only species 
from S) a smooth function 

R> x R| -> R s 
(K k ,x) H> Rl k {x) 

such that 

• for K k > 0, the function R^ k (x) is a rate function for the reaction y k — > y' k that 



satisfies Assumption 2.6, and 



when K k = 0, then R^ k (x) is the zero function. 
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2. Let (5 be a chemical reaction network, and let K, be a "parametrized family of kinetics on 
\S\ species. Then & is said to admit multiple tC steady states or is /C-multistationary if 
there exist kinetics {R^, 1 *} arising from K, and a stoichiometric compatibility class V such 
that the resulting system Q has two or more positive steady states in V . Moreover, such 
a network is said to admit bistability if such steady states can be found that are stable. 

As noted above, an important family of kinetics fC is that of mass-action kinetics; in this 
case, (3 admits multiple mass-action steady states if there exist rate constants Kk £ M>o and a 
stoichiometric compatibility class V such that the mass-action system Q admits at least two 
positive steady states in V . 

Example 2.10. We again consider the CFSTR Q examined in Examples 2.5 and 2.8. Recall 
that for a CFSTR, the unique stoichiometric compatibility class is the nonnegative orthant: 
here, V = M> . Therefore, our CFSTR ([I]) admits multiple positive mass-action steady states 
if and only if there exist reaction rate constants K\, k 2 , ■ • • > K io £ ^>o such that the differential 
equations (fsl) have at least two positive steady states. Indeed, the CRN Toolbox flOf determines 
that when the mass-action system takes the following rate constants: 

(Kl, K 2 , Ki, K 5 , K 6 , K 7 , Kg, Kg, Kl ) = 

(1, 1, 1, 1, 41774.858, 1, 2.5081 * 10" 4 , 7.3335 * 1CT 3 , 1.1614 * 1(T 4 , 7.5610 * 10" 5 ) , 

there are two steady states: 

x* = (63.143335,136.35902,41577.356) and 
x** = (25473.839, 1007.5644, 15295.454) . 

2.3 Classification of multistationary one-reaction CFSTRs 

We now recall the following theorem, due to the first author: 

Theorem 2.11 ([IS]). 1. Consider a CFSTR which contains only one non-flow reaction: 

aiXi + a 2 X 2 + ■■■ + a s X s 61X1 + b 2 X 2 + ■■■ + b s X s , 

where a^,bi > 0. Then the CFSTR admits multiple positive mass-action steady states if 
and only if b i>a .ai > 1- Moreover, these multistationary CFSTRs admit nondegen- 
erate steady states. 

2. Consider a CFSTR in which the only non-flow reactions consist of a pair of reversible 
reactions: 

aiXi + a 2 X 2 + ■■■ + a s X s *=> b x Xx + b 2 X 2 + ■■■ + b s X s , 

where ai,bi > 0. The CFSTR admits multiple positive mass-action steady states if and 
only if the following holds: 

a,[ > 1 or h > 1 . 

i: bi>ai i: a,i>bi 

Moreover, these multistationary CFSTRs admit nondegenerate steady states. 
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The current work was motivated by the question of whether a similar theorem exists for 
the class of CFSTRs that consists of networks with two reversible nonflow reactions and their 
sub-CFSTRs. 



3 Lifting multistationarity from subnetworks 

Consider the following question: if a subnetwork N of a network G admits multiple positive 



steady states, then does G as well? Theorem 3.1 asserts that the answer to this question is 



'yes,' provided that the steady states are nondegenerate and the two networks share the same 
stoichiometric subspace (note that the stoichiometric subspace of N is always contained in that 
of G). The proof lifts each steady state x* of N to a nearby steady state of G. 

Theorem 3.1. Let N be a subnetwork of a chemical reaction network G such that they have 
the same stoichiometric subspace: Sn = Sq- Let K, be a parametrized family of kinetics on the 
species of G. Then the following holds: 

• If N admits multiple nondegenerate positive K, steady states, then G does as well. Addi- 
tionally, if N admits finitely many such steady states, then G admits at least as many. 

• Moreover, if N admits multiple positive exponentially stable steady states, then G does as 
well. Additionally, if N admits finitely many such steady states, then G admits at least 
as many. 

We note that our theorem is similar to Theorem 2 in work of Craciun and Feinberg [4j; 
their theorem allows multiple steady states to be lifted from an 'entrapped species' network 
(that is, only certain species are in the outflow) to the corresponding 'fully diffusive' network 
(all species are in the outflow). In addition, their theorem is stated as a contrapositive version 
of ours. Our proof of Theorem |3. 1| makes use of the following homotopy theory result, which is 
a modified form of Theorem 1.1 in Craciun, Helton, and Williams [7]. 

Lemma 3.2. Let S C K" be a vector subspace, let V C R n be a polyhedron contained in an 
affine translation of S, and let Q C int("P) be a bounded domain in the relative interior of V . 
Assume that g\ : $7 — >• S , for A G [0, 1], is a continuously-varying family of smooth functions 
such that 

1. for all A G [0, 1], g\ has no zeroes on the boundary of VL, and 

2. for A = and A = 1, Im.dg\(x) = S for all 

Then the number of zeroes of go in f2 equals the number of zeroes of gi in Q. 
We now prove Theorem |3.1| 



Proof of Theorem 3.1 First, note that the network G and its subnetwork iV must have the 
same set of species S in order for their stoichiometric subspaces to coincide. We let S denote 
the shared stoichiometric subspace: S := Sg = Sjy- Now, let 1Z' denote the set of reactions of G 
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that are not in N: TZq = TZ-n U TZ'. We now assume that the subnetwork N admits multiple 
nondegenerate positive steady states; that is, there exist rate constants k\, k* 2 , . . . , kL, G M>o 
such that there exist distinct, stoichiometrically compatible, nondegenerate positive steady 
states x* and x** of the chemical reaction system (N,n*) arising from fC. Write /jv for the 
differential equations of (N,k*). Now x* is a nondegenerate steady state, so there exists a 
relatively open ball f2 around x* in the interior of V such that (1) x* is the unique steady state 
(zero of /tv) in O, and (2) fcndfjs?(x) = S for all x G f2. Note that (2) can be accomplished 
because the non- vanishing of a determinant is an open condition and because the matrix df^ix) 
varies continuously in x. 

For any vector of reaction parameters n G M^q , we define the following the following family 
of functions for < A < 1: 

gl{x) := Mx ) + ^(y' k -y k )R^( x ) . 
hew 

It follows that g\(x) gives the differential equations Q of the chemical reaction system arising 
from the network G and the following reaction parameters with respect to the kinetics tC: 

(«*,«) := ( K Kl,Kl,...,K* nNl ,XK 1 ,\K 2 ,...,X^ n >^ £M l ^ al . (9) 

Note that Qq (x) = /n(x). Next, by continuity in k and the compactness of the boundary of Q, 

there exists a vector of reaction parameters G ' such that for all < A, the function g*f (x) 
has no zeroes on the boundary of f2. By continuity in A, and by scaling smaller if nec essary 

allows 



3.2 



we may assume additionally that Im d x g*f (x) = S for all x G il. Therefore, Lemma 
us to conclude that the chemical reaction system (G,(k*,Xk')) has a nondegenerate steady 
state in the ball Q, for all < A < 1. 

We now complete the proof by repeating the argument with 5;** , taking care that the ball 
around x** does not intersect that of x*; we replace by a scaled-down version (/jlk^ for some 
< n < 1) if necessary. It follows that (G, (k*, Ak')) has at least two nondegenerate steady 
states. The case of three or more nondegenerate steady states generalizes in a straightforward 
way. For the stability result, we simply note that the eigenvalues of a matrix vary continuously 
under continuous perturbations (in this case, arising from the parameter A). □ 



Remark 3.3. One application of Theorem 3.1 is that it provides the mathematical justifica 



tion for the technique of Siegal-Gaskins et al. which establishes bistability in the mass- action 
setting by way of 'network ancestry' J23^ . In their examination of 40, 680 small gene regulatory 
networks, 14, 721 initially were established to be bistable by the implementation of Advanced De- 
ficiency Theory in the CRN Toolbox JiOj/ . and an additional 22, 050 were classified as bistable by 
virtue of containing one of the 14, 721 bistable networks as a subnetwork ('ancestor') such that 
both networks have the same stoichiometric subspace. A similar approach is taken by Conradi et 
al. for lifting multiple steady states from certain subnetworks called 'elementary flux modes ' J1J/. 
We note that their criterion for lifting steady states does not require that the stoichiometric 
subspaces of the network and its subnetwork to coincide JH Supporting Information]. 



10 



The next example illustrates why the hypothesis of nondegeneracy is required in Theo- 



rem 



3.1 A larger such example appears in the work of Craciun and Feinberg [H §6]. 



Example 3.4. Consider the following (non-CFSTR) network: 




A + C 



2B 



(10) 



The CRN Toolbox ]10$ determines that network (10) does not admit multiple positive mass- 
action steady states, but the following subnetwork does admit multiple degenerate positive steady 
states: 



A^B^ 



C 



A + C ^ 2B . 



(11) 



In fact, it is straightforward to verify that steady states exist for network (11) if and only if 
k-i = k>2, and in this case, each two-dimensional compatibility class contains an infinite one- 
dimensional set of degenerate steady states. 

One way for a network and its subnetwork to share the same stoichiometric subspace is 
for the subnetwork to be obtained by making some reversible reactions irreversible. Thus, 



Theorem 3.1 yields the following corollary. 



Corollary 3.5. For a chemical reaction network N , let )C be a parametrized family of kinetics 
on the species of N. Let G be a network obtained from N by making some irreversible reactions 
of N reversible. Then if N admits multiple nondegenerate positive K, steady states, then G does 
as well. 



The next corollary states that Theorem 3.1 allows multiple positive steady states to be 
lifted from a sub-CFSTR to a fully open CFSTR. Therefore, the set of minimal multistationary 
CFSTRs (with respect to the subnetwork relation) completely defines the set of all fully open 
multistationary CFSTRs: a fully open CFSTR admits multiple steady states if and only if it 
contains as a subnetwork one of these minimal CFSTRs. This result will be useful in our 
classification of small multistationary CFSTRs in Section [7j 

Corollary 3.6. Let N be a sub-CFSTR of a fully open CFSTR G, and let K, be a parametrized 
family of kinetics on the species of G. Then, if A admits multiple nondegenerate positive K. 
steady states, then G does as well. 



Proof. Assume that the species of A are X\, A2, . . . , X Sl and the species of G are X%, X2, ■ 



X 



S1+S2- 



Let N' be the CFSTR obtained from N by appending the flow reactions ^ A 



Sl+l) 



^ A Sl+ 2, . . . , ^» X Sl J rS2 for all species of G that are not in N. Clearly, N' is a subnetwork of 
G, and they share the same stoichiometric subspace, namely, M Sl+S2 . By applying Theorem 3.1 
to A' and G, we see that if A' admits multiple nondegenerate positive steady states, then G 
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does as well. Therefore, it remains only to show that ./V admits multiple nondegenerate positive 
steady states if and only if N' does. 

Consider any outflow rate parameter K°^\^ > for one of the new outflow reactions. Then 

by Assumption 

x Sl+ i from 



2.6 



the rate function R X ^ t + '^ v (x) depends only on x Sl+ i and is increasing in 

K s 1 +i 

— R^out > • • • j ^si+i— 1 1 0, , . . . X Sl -\- S2 ) • 

si+i 



2.6 



implies that R°^ Xai+l (x) is a 

K s 1 +i 



As for the corresponding inflow rate function, Assumption 
positive constant function, and this constant depends only on the parameter K™ +i and is in fact 

increasing in this parameter for sufficiently small values, with R 81+1 (x) = 0. Thus, we can 
choose a sufficiently small inflow parameter > such that there exists a positive value 

x* Sl+ i > for which the rate functions are equal at x when x Sl+ i = x* i+i . 

Therefore, it follows that x* = (a;*, x\, . . . ,x* Sl ) £ M> is a nondegenerate positive steady 
state of the system (iV, («i, K2, • • • , k\k n \)) if and only if the concentration vector 

^aj 1; x 2 , • • • , x Sl , x Sl+ i, . . . , x Sl+S2 j t i& >0 
is a nondegenerate positive steady state of the system 

(N', (ki,K2, • • • , K\K N \> K n+1' K si+1 ■ ■ ■ J K °H-s 2 )) ' 

where the rates and fc™^ are chosen as described above. This completes the proof. □ 

4 Lifting mass-action multistationarity from embedded CFSTRs 



Corollary 3.6 stated that multistationarity can be lifted from sub-CFSTRs; in this section, we 
generalize the result to the case of embedded CFSTRs in the mass-action setting (Corollary |4.6[ ). 

We first need the following generalization of inflow/outflow reactions in order to allow for 
reactions such as A £5 2A which also have a mass-action steady state at x~X = 1 when the two 
reaction rate constants are equal. 

Definition 4.1. A mass-action flow- type subnetwork for a species Xi of a chemical reaction 
network G is a nonempty subnetwork N of G such that 

1. the reactions in N involve only species Xi, and 

2. there exists a choice of reaction rate constants k* for the reactions r S IZm of N such 
that for the resulting mass-action system of this subnetwork N , x~l = 1 is a nondegenerate 
steady state. 



The following theorem is analogous to Theorem 3.1 



Theorem 4.2. Let N be an embedded network of a network G such that 
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1. the stoichiometric subspace of N is full- dimensional: Sjy = 1$) Sn \, and 

2. for each species Xi that is in G but not in N, there exists a mass-action flow-type sub- 
network of G for Xi . 

Then the following holds: 



If N admits multiple nondegenerate positive mass-action steady states, then G does as 
well. Additionally, if N admits finitely many such steady states, then G admits at least 
as many. 

Moreover, if N admits multiple positive exponentially stable mass-action steady states, 
then G does as well. Additionally, if N admits finitely many such steady states, then G 
admits at least as many. 



The proof of Theorem 4.2 requires the following lemma, which states that for certain simple 
embedded networks obtained by removing only one species, each nondegenerate steady state u 
can be lifted to a steady state of the larger network that is near (u, 1). 

Lemma 4.3. Let G be a chemical reaction network with s species denoted by X\,X 2 , ■ ■ ■ ,X S , 
and let N be an embedded network of G with s — 1 species X\,X 2 , . . . ,X s _i such that N is 
full- dimensional: Sn = M s_1 . Assume that the reactions of G and the reactions of N can be 
written as, respectively, TZq = {Ri, R2, ■ ■ ■ , Rm, Rm+i, ■ ■ ■ , Rm+n} and TZn = {Ri, R2, ■ ■ ■ , Rm} 
such that: 

1. for i = 1,2, ... ,m, the reaction Ri of N is obtained from the corresponding reaction Ri 
of G by removing species X s , and 

2. all remaining reactions of G, namely {R m+ \, R m+2 , ■ ■ ■ , R m +n}, together form a mass- 
action flow-type subnetwork for the species X s . 

For a choice of rate constants k\, k\, . . . , K* m > 0, let £ (N, {k*, k^, . . . , Km}) denote a finite set 
of nondegenerate positive mass-action steady states of the system arising from N and the k*. 
Then for sufficiently small e > 0, there exist reaction rate constants K m+1 , K m+2 , ■ ■ ■ , K m +n f or 
the flow-type subnetwork ofG such that for allu € X (N, {n^, k\, . . . , K m }), there exists a nonde- 
generate positive mass-action steady state u of the system arising from G and k\, k^, ■ ■ ■ , K m + n 
with \u — (u, 1)| < e. Additionally, if u is exponentially stable, then u is as well. 

Proof. Fix a choice of reaction rate constants k*, k\, . . . , n m , and let S := S (TV, {k\, ■ ■ ■ , K m }) 
be as in the statement of the lemma. 

We view G = N(jM as the disjoint union of two subnetworks, one which consists of the reac- 
tions R±, i?2j ■ • • > Rm, which we denote by N, and the second which consists of R m +±, R m +2, ■ ■ ■ > 
R m+n , which we denote by M. As M is a mass- action flow- type subnetwork for the species X s , 
there exist rate constants K m+ i, K m+ 2, • • • , n m +n > such that the resulting mass-action ODE 
system, denoted by /m(^s) 5 has a nondegenerate steady state at = 1. 
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Next, we denote by f^{x) the mass-action ODE system ([7]) arising from the subnetwork N 
and the fixed rate constants k*, k^, . . . , K* m . Consider the following map from M>o x R> to M s : 

f G (k,x) := (fff tl (x), f^ y2 {x), ... Jft^ix), ffi s (x) + kf M (xs)) • (12) 

(Note that f^. denotes the i-th coordinate function of f^.) It follows that fc(k,x) denotes 
the mass-action ODEs for the network G with respect to the rate constants 

K-y, . . . , K m , /cK m _|_l, A)K m -|_2, . . . kK m + n . 

We scale the last coordinate of fc(k,x) by 1/k and make the substitution 6 = 1/k to obtain: 

F G (5,x):=(f Rl (x), ... f M (x s )) + (0,...,0,<S/# j ») (13) 

=:h(x) + s{0,...,0,f Rs (x)^ , 

where h(x) is defined by the second equality. Hence, it suffices to prove that for sufficiently 
small e > and for all u G E, there exists a 5 > such that there exists a nondegenerate zero 
u of Fg(5, x) with \u — (u, 1)| < e. 

Fix ueS. We now claim that h has a nondegenerate zero at (u, 1). The final coordinate 
of /i satisfies h s (u, 1) = /m(1) = by construction. As for the remaining coordinates z = 
1, 2, . . . , s — 1, we compute 

hi(u,l) = = /*(«) = . (14) 



We now explain the second equality in (14). When the reaction i?-,- in iV is yj — > y'j then 
the reaction Rj, given by yj — > y'j, is such that the projection of y~j onto the first s — 1 

coordinates is yj and similarly for y^.. Thus, the reaction vector y'- — y~j projects to y'j — yj and 
U2, . . . , u s _i, l) Vj = (u\, U2, ■ ■ ■ , u s -i) Vj '. Finally, (u, 1) is nondegenerate, because dh(u,l) 
is an s X s-matrix in which the upper-left (s — 1) x (s — l)-submatrix is the nonsingular matrix 
df N (u) and the bottom row is (0, 0, ... , 0, £-(1)) with f^(l) + by hypothesis. 

As (u, 1) is nondegenerate, there exists a constant e(u) > such that the resulting e(u)- 
neighborhood of (u, 1), which we denote by is such that (1) f2 is in the positive orthant M^q, 
(2) (u, 1) is the unique zero of h in fi, and (3) dh{x) is nonsingular for all x 6 f2. Consider 



again the function Fg(S, x) defined in (13), and note that F(0,x) = h(x). By continuity in 5 
and the compactness of the boundary of Q, there exists 5(u) > such that for all < 5 < 
6(u), the function Fg(5,x) has no zeroes on the boundary of O. Again by continuity and by 
decreasing S(u) if necessary, we may assume that dF{5, x) (the matrix of partial derivatives with 
respect to the x\, X2, • ■ • , x s ) is nonsingular for all < S < 5(u) and for all x G 17. Therefore, 
Lemma |3.2| allows us to conclude that Fq(6,x) has a unique nondegenerate zero u in Q (that 
is, \u - (u, 1)| < e(u)) for all < 5 < 5(u). 

Now let e* be the minimum of all such e(u), where u G S. Additionally, we decrease e* if 
necessary so that the resulting e-neighbor hoods of the points u do not intersect. The lemma 
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now follows with the e* as a cut-off: given any < e < e*, the above arguments for each u can 
be made using e in place of e(u). Taking the minimum, denoted by 5* , of the resulting cut-offs 
5(u), we obtain nondegenerate zero u of Fg(S* , x) such that \u — (u, 1)| < e. 

For the stability result, the eigenvalues of a matrix vary continuously under continuous 
perturbations (in this case, arising from the parameter 5). □ 



Remark 4.4. In the proof of Lemma 4-3, the s- dimensional dynamical system (12) may be 
represented by 

x = In(. x ) + k 9M{x s ) , 

where guii^s) '■= (0, . . . , 0, /m(^s))- When k is sufficiently large and x s is in the domain of 
attraction of ~x~2 = 1 with \gM{x s )\ ~ 0(1), then x ~ kgM{x s ), which has dynamics close 
to the one- dimensional system x s = /c/a/(x s ). However, when x s is close to Xs = 1 with 
\gM(x s )\ ~ 0(l/k 2 ), then x ~ fft(x), the dynamics of which are effectively those of an (s — 1)- 
dimensional system. Thus by choosing k large enough, we achieve a time-scale separation: on 
the fast time-scale, the dynamics are close to a one- dimensional system and on the slow time- 
scale, the dynamics are close to an (s — 1) -dimensional system. Thus, we can lift the steady 
states from the smaller system to the full system. 



We can now prove Theorem 4.2 



Proof of Theorem \4-S\ We begin by reducing to the case that G has only one species that 
iV does not have: if G has more than one additional species, we can lift multistationarity 
'one species at a time.' Now denote the species of N by X\, X2, ■ ■ ■ , X s -i and the species 
of G by X 1 ,X 2 ,...,X S . Denote the reactions of N by y 1 -)■ y[,y 2 -4 y' 2 ,...,y m -)• y' m , 
where yi,y[ S <^>o • As N is an embedded network of G, we can write the reactions of G as 

T^-G = {Rl> R2, ■ • • j Rmi Rm+l, Rm+2, ■ ■ ■ i Rm+m Rm+n+li ■ ■ ■ i Rm+n+p} Such that: 

1. for % = 1,2, ... ,m, the reaction Ri of N is obtained from the corresponding reaction Ri 
of G by removing species X s , and 

2. the reactions in {R m+ i, R m +2, ■ ■ ■ , Rm+n} form a mass-action flow-type subnetwork for 
the species X s . 

We now let G' denote the subnetwork of G that consists of the reactions: R\, R 2 , ■ ■ ■ , R m , 



Rm+i, Rm+2, ■ ■ ■ , Rm+n- Lemma 4.3 applies to this network G' and its embedded network N, 
so G' admits at least as many nondegenerate positive mass-action steady states as N (and 
similarly for exponentially stable steady states). Next, Q is a subnetwork of G that shares 



the same stoichiometric subspace (namely, M s ), so by Theorem 3.1, G admits at least as many 
nondegenerate positive mass-action steady states as G' (and similarly for exponentially stable 
ones), so this completes the proof. □ 



We now illustrate the necessity of the hypothesis 2 of Theorem 4.2 
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Example 4.5. Consider the following (non-CFSTR) network G, which is adapted from a sim- 
ilar network that appears in work of Feinberg }12$ : 



+± A 3At±2A + B 

A straightforward calculation reveals that G has a unique mass-action steady state, namely 
(xa, xb) = ^^f)- I n f ac ^! despite the fact that A participates in a non-flow reaction, the 
steady state value of xa is the same as it would be when considering only the flow subnetwork 
^ A. Now consider the following embedded network N obtained by removing the species B: 

±5 A 3A ±3 2A 



We see that N satisfies the conditions of Theorem \2.11, so N admits multiple mass-action 



steady states. Note that N is an embedded network of G, but its multiple steady states can not 



be lifted to G; Theorem 4-2 does not apply because G does not contain a flow-type subnetwork 
for the species B. 

On the other hand, N is an embedded network of the following network G' : 

B^O^i 3A^2A + B 



which does contain a flow-type subnetwork for the species B. So, Theorem 4-2 does apply and 
thus we conclude that G' admits multiple steady states. 

We now have an analogue of Corollary |3.6| 



Corollary 4.6. Let N be an embedded CFSTR of a fully open CFSTR G. Then, if N admits 
multiple nondegenerate positive mass-action steady states, then G does as well. 



Proof. This follow directly from Theorem 4.2 after noting that hypothesis 2 of the theorem is 



satisfied by the inflow/outflow reactions ^ Jj. □ 

5 CFSTR atoms of multistationarity 

In the previous section, we saw that a CFSTR is multistationary in the mass-action setting if and 
only if an embedded CFSTR is multistationary; now we call the minimal such networks 'atoms of 
multistationarity. ' In Section [7j we will classify certain two-reaction atoms of multistationarity 



(see Corollary 7.2). 



Definition 5.1. 1. A fully open CFSTR is a CFSTR atom of multistationarity if it admits 
multiple nondegenerate positive mass-action steady states and it is minimal with respect 
to the embedded network relation among all such fully open CFSTRs. 

2. A fully open CFSTR G is said to possess a CFSTR atom of multistationarity if there 
exists an embedded network N of G that is a CFSTR atom. 
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We now restate Corollary 4.6 in the following way, which motivates the above definition and 
suggests that compiling a list of atoms is desirable. 



Corollary 5.2. A fully open CFSTR possesses a CFSTR atom of multistationarity if and only 
if it admits multiple nondegenerate positive mass-action steady states. 

Proof. The reverse direction is clear: a multistationary CFSTR is either itself a CFSTR atom 



of multistationarity or contains one. The forward direction is Corollary 4.6 □ 



We also can rephrase Theorem 2.11 in the following way: 



Corollary 5.3. A one-reaction CFSTR is a CFSTR atom of multistationarity if and only if it 
consists of one non-flow reaction and that non-flow reaction has one of the following two forms: 

ai X a 2 X , or X + Y -> b 1 X + b 2 Y , (15) 

where a 2 > ai > 1, or, respectively, b\ > 1 and b 2 > 1. A one-reaction CFSTR possesses one 
such CFSTR atom of multistationarity if and only if it admits multiple nondegenerate positive 
mass-action steady states. 

We end this section by posing the following questions: 

1. Is there a good characterization of CFSTR atoms of multistationarity? For instance, even 
though there are countably infinitely many one-reaction CFSTR atoms, Corollary 1 5 . 3| gives 
a simple characterization of all such one-reaction atoms. In particular, a one-reaction atom 
contains at most two species, and furthermore each of these atom types is characterized 



by exactly two parameters, (01,02) or (61,62) in equation (15). 



2. Is there a good notion of 'atom of multistationarity' outside of the CFSTR setting? If 
so, then a CFSTR atom might contain as an embedded network, a more general atom, 
which is obtained by removing some flow reactions and possibly more reactions. For 
example, we can remove the outflow reaction A — > from the CFSTR atom arising 
from A — > 2A A + B — > (see the top of Figure [3] in the next section) and maintain 
multistationarity, but removing B — > destroys multistationarity. 

Beginning in the next section, we will give a partial answer to the first question above for 
two-reaction CFSTRs. 

6 Enumeration of reversible bimolecular two-reaction CFSTRs 

The remainder of this work is dedicated to answering the following question: 

Question 6.1. Which bimolecular two-reaction fully open CFSTRs admit multiple positive 
mass-action steady states? 
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By bimolecular we mean that each complex contains at most two molecules: the complexes 0, 
A, 2A, and A+B are permitted, but 2A+B is not. A two-reaction CFSTR refers to a CFSTR in 
which the non-flow reactions consist of two pairs of reversible reactions, one reversible reaction 
and one irreversible reaction, or two irreversible reactions. For instance, the three CFSTRs Q, 
Q, and Q in Example 2.5 are among the bimolecular two- reaction CFSTRs for which we 



would like to answer Question |6.1| Let us note that reactions of the form ^ 2 A or ^5 A + B 
(or any of the directed versions) are considered non-flow reactions. Finally, if we define two 
networks to be equivalent if there exists a relabeling of the species that transforms the first 
network into the second network, we aim to list only one network from each such equivalence 
class. For example, the two CFSTRs in which the non-flow reactions are 2A ^— A + B <— A 
and C — > B + C — > 2C, respectively, are both in the same equivalence class. 

Note that it is sufficient to enumerate the possible non-flow subnetworks of our CFSTRs of 
interest; for example, if the non-flow subnetwork is 

2A^A + B^A, (16) 

then the corresponding CFSTR is obtained by including the flow reactions for species A and B. 



In addition, Corollary 3.5 implies that a non-reversible CFSTR (for example, the one arising 
from (16)) does not admit multiple nondegenerate positive steady states if the corresponding 
reversible CFSTR (for example, the one arising from 2A ^ A + B ^ A) does not. Therefore, 
we will proceed to answer Question |6.1| by completing the following steps: 



1. Enumerate all reversible bimolecular two-reaction networks. 

2. Determine which of the fully open CFSTRs arising from networks in Step 1 admit multiple 
positive mass-action steady states. 

3. Of those reversible CFSTRs that admit multiple positive steady states which were found 
in Step 2, determine which sub-CFSTRs admit multiple positive steady states. 



The current section describes how we performed Step 1 (see Algorithm 6.4), and in Section [7| 
we explain how we completed Steps 2 and 3. 

6.1 The total molecularity partition of a chemical reaction network 

We now explain how a network defines a 'total molecularity partition'; two-reaction networks 



will be enumerated by these partitions in Algorithm 6.4 Recall that a partition of a positive 
integer m is an unordered collection of positive integers that sum to m; by convention, we write 
the partition as (mi, m^, • • • , m n ), where the parts nii are weakly decreasing: m\ > mi > ■ ■ ■ > 
m n > 1. Partitions of m = 4, 5, 6, 7, 8 are listed (partially) in Table [T] 

Example 6.2. Let us rewrite the network 2A ^A + B^A + C as two separate reversible 
reactions: 

2A^A + B A + B^A + C. (17) 
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Table 1: Here we (partially) list the partitions of m = 4, 5, 6, 7, 8 in lexicographic order. The numbers 
of partitions, which are known as the Bell numbers, are listed in the last column. 



m 


Partitions of m 


# of partitions 


4 


(4), (3, 1), (2, 2), (2, 1, 1), (1, 1, 1, 1) 


5 


5 


(5), (4, 1), (3, 2), (3, 1, 1), (2, 2, 1), (2, 1, 1, 1), (1, 1, 1, 1 , 1) 


7 


6 


(6), (5, 1), (4, 2), (4, 1, 1), (3, 3), (3, 2, 1), ... , (1, . . . , 1) 


11 


7 


(7), (6, 1), (5, 2), (5, 1, 1), (4, 3), (4, 2, 1), . . . , (1, . . . , 1) 


15 


8 


(8), (7, 1), (6, 2), (6, 1, 1), (5, 3), ... , (1, 1, 1, 1, 1, 1, 1, 1) 


22 






60 



Counting the number of times each species appears (where we take into consideration the sto- 
ichiometric coefficients), we see that species A appears 5 = 2 + 1 + 1 + 1 times, B appears 
2 = 1 + 1 times, and C appears 1 time. Definition 6.3 will say that the 'total molecularities' of 
species A, B, and C are, respectively, 5, 2, and 1. In addition, the 'total molecularity partition' 
of network ( 17) will be (5, 2, 1), which is a partition of the integer 8 = 5 + 2 + 1. Similarly, the 
total molecularity partition of the network 2 A t> A + B ^ A is (5, 2), a partition of 7. 

The definition of total molecularity first appeared in |17j . 

Definition 6.3. 1. For a reversible network, the total molecularity of species Xj refers to 
the sum over all pairs of reversible reactions of the sum of the stoichiometric coefficients 
of Xj in the reactant and in the product: 

TM(Xj) := a i +b ^ 
El=i ai^ELi hXi e n> 

where 1Z! denotes all pairs of reversible reactions Y^t=i a iXi ^ biX{. 

2. For a reversible network, the total molecularity partition is the partition defined by the 
multiset of total molecularities of all species: 

{TM(X0, TM(X 2 ), ... , TM(X| 5 |)} . 

Note that for reversible bimolecular two-reaction networks, the total molecularity partition 
is of an integer m £ {4, 5, 6, 7, 8}. 



6.2 Algorithm for enumerating networks 

We now present the algorithm we used for enumerating reversible bimolecular two-reaction 
networks. 

Algorithm 6.4 (Algorithm for enumerating reversible bimolecular two-reaction networks). 
Step One. List partitions of m = 4, 5, 6, 7, 8. 
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Step Two. For each partition (mi, TO2, ■ ■ ■ , m n ), list (with repeats) all reversible bimolecu- 
lar two-reaction networks in which species X\ has total molecularity m\, species X2 has total 
molecularity mi, and so on. 

Step Three. Remove networks that contain trivial reactions, networks that contain repeated 
reactions, and decoupled networks. 

Step Four. Remove redundant networks: keep exactly one representative from each equiva- 
lence class of networks. (Recall that two networks are equivalent if there exists a relabeling of 
the species that transforms the first network into the second network.) 



As we see in Table [2j Algorithm 6.4 yields 386 reversible bimolecular two-reaction networks. 
In Section [7| we determine which of the 386 CFSTRs admit multiple positive steady states. 



Let us now elaborate on our implementations of Steps Two through Four of Algorithm 6.4 



In order to list all reversible bimolecular two-reaction networks that have a given partition 
(mi, m2, . . . , m n ) (Step Two), we made use of a psuedo-species Xq. Namely, any network with 
partition (m±,m2, • • • , m n ) arises from placing (m — m\ — 1112 — ■ ■ ■ — m n ) copies of species Xq, 
mi copies of X\, mi copies of X2, and so on in the eight boxes in the following diagram: 



□ + □ 



□ + □ 



□ + □ 



□ + □ 



For example, A + A ±5 A + B A + B ±5 X + A defines the network 2 A £5 A + B 



A. Clearly, this procedure will yield all networks, but certain trivial networks (such as one 
with repeated reactions) will appear, and additionally each network will appear more than 
once. Accordingly, trivial networks are removed in Step Three of Algorithm |6.4[ and Step Four 
keeps only one representative from each equivalence class of networks. Step Four is the most 
computationally expensive part of our enumeration. For each network remaining at the end of 
Step Three, we generated the equivalence class of networks obtained by performing a relabeling 
of the species. Two networks are equivalent if and only if they generate identical equivalence 
classes of networks. We removed extra copies of equivalent networks at the end of Step Four. 



6.3 The enumeration of small networks of Deckard, Bergmann, and Sauro 

A related (and much larger: over 47 million) enumeration of small bimolecular networks was 
undertaken by Deckard, Bergmann, and Sauro [Sj. Their work enumerated small networks by 
the number of directed reactions and by the number of species. So, the network 2A <— A+B <— A 
falls in their list of networks containing two directed reactions and two species, and the network 
2A ^ A + B -<r- A + C is & network containing three directed reactions and three species. Also, 
their enumeration did not include seemingly unrealistic chemical reactions involving the zero 
complex (such as — > 2A or A + B) or reactions in which some species appears in both 
the reactant complex and product complex of a reaction (such as ^4— > A + B or A — > 2A). We 
remark that from this enumeration of networks by Deckard, Bergmann, and Sauro, the work of 
Pantea and Craciun sampled networks to compute the fraction that pass the Jacobian Criterion 
[201 Figure 1]. 



20 



7 Classification of multistationary two-reaction CFSTRs 



The main result of this section is the following theorem, which completely answers Question |6.1| 

Theorem 7.1. Of the 386 reversible, bimolecular, two-reaction fully open CFSTRs, exactly 
35 admit multiple positive mass-action steady states. Moreover, each of these 35 networks 
admits multiple nondegenerate positive steady states. Furthermore, each such network contains 
a unique minimal multistationary subnetwork. The poset (partially ordered set) of these 35 
directed subnetworks, with respect to the embedded network relation, has 11 minimal elements, 
which are the bimolecular two-reaction CFSTR atoms of multistationarity. 



An immediate corollary of Theorem 7.1 is the following: 



Corollary 7.2. 1. A bimolecular, two-reaction fully open CFSTR admits multiple nonde- 
generate positive mass-action steady states if and only if it contains as a sub-CFSTR one 
of the 35 minimal such subnetworks, which are displayed in Figure [J| 

2. A bimolecular, two-reaction fully open CFSTR admits multiple nondegenerate positive 
mass-action steady states if and only if it contains as an embedded network one of the 11 
CFSTR atoms which are marked in bold/red in Figure^ 

3. If a fully open CFSTR (not necessarily bimolecular and having any number of reactions) 
G contains one of the 35 minimal CFSTRs mentioned above as a sub-CFSTR or contains 
one of the 11 atoms as an embedded network, then G admits multiple nondegenerate 
positive mass-action steady states. 



Note that part 3 of Corollary |7.2| makes use of Corollaries 3.6 and 4.6 



Example 7.3. Among the 35 reversible CFSTRs in Theorem 7.1 that admit multiple steady 
states, one is the network ([I]) which we first saw in Example 2.5- it arises from the network 
2 A ^-A + B^A + C. The unique minimal multistationary sub-CFSTR is the directed sub- 
network ([2]) obtained by removing two reactions: 2A <— A + B <— A + C . Finally, there is 
a multistationary embedded CFSTR ^ obtained by removing species C, namely, the CFSTR 
arising from 2 A ^— A + B A, that is one of the 11 atoms. In other words, no further embed- 
ded CFSTR is multistationarity. The directed subnetwork ^ and the embedded network Q 
appear in the lower left of Figure [3| 



Sections 7.1 through |7.3| provide the proof for Theorem 7.1 



7.1 Ruling out multistationarity by the Jacobian Criterion 

Recall from [31 HI [6] that the Jacobian Criterion is a method for ruling out multistationarity. 
A CFSTR is said to pass the Jacobian Criterion if all terms in the determinant expansion of 
the Jacobian matrix of its mass-action differential equations ([7]) have the same sign. Craciun 
and Feinberg proved that if a CFSTR passes the Jacobian Criterion, then it does not admit 
multiple positive steady states. In earlier work, the current authors proved that if the total 
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# 


# 


# 




reversible bimolecular two-reaction networks 


#of 


with 


that fail 


with 


m 
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networks 


TM > 2 


Jac. Crit. 


MSS 


4 


(0,2,2,5,3) 


12 


2 








5 


(1,4,7,8,10,9,2) 


41 


20 


8 


1 


6 


(0,3,6,9,7,23,12,9,23,12,3) 


107 


60 


31 


5 


7 


(0,1,3,4,5,13,7,9,13,26,8,12,15,7,1) 


124 


89 


55 


15 


8 


(0,0,0,1,1,3,2,1,5,4,9,4,7,8,13,12,3,5,11,9,3,1) 


102 


73 


48 


14 






386 


244 


142 


35 



Table 2: Here we list the number of reversible bimolecular two-reaction CFSTRs by partition. The 
order of partitions is the lexicographic order (as in Table [lj . In bold are the 142 networks for which 
multistationarity is ruled out because each part of the corresponding partitions (the total molecularity of 
a species, denoted by "TM") is no more than two [17]. Of the remaining 244 networks, an additional 102 
networks pass the Jacobian Criterion (those with total molecularity at most two also pass the Jacobian 
Criterion). For the remaining 142 networks, the CRN Toolbox [TU] determined that precisely 35 admit 
multiple positive mass-action steady states ("MSS"). 



molecularities of all species are at most two, then the CFSTR passes the Jacobian Criterion [17] . 
Accordingly, any two-reaction networks that arise from the 19 partitions (of 4, 5, 6, 7, or 8) in 
which all parts are at most two automatically pass the Jacobian Criterion; these 142 networks 
are marked in bold in Table [2j Of the remaining 244 = 386 — 142 networks, an additional 102 
networks pass the Jacobian Criterion. 

7.2 Applying the CRN Toolbox to classify reversible two-reaction networks 

For the remaining 142 reversible networks that do not pass the Jacobian Criterion, we applied 
the CRN Toolbox |10j . This was performed in an automated fashion by using Autolt code p] 
provided by Dan Siegal-Gaskins. We find that exactly 35 admit multiple positive mass-action 
steady states and the remaining 107 do not. For each of the 35 multistationary CFSTRs, the 
Toolbox gave an instance of rate constants, two positive steady state values, and the correspond- 
ing eigenvalues. In all cases but one, the nondegeneracy of these steady states was evident from 
the eigenvalues. In the remaining case, in which one steady state was degenerate, we found 'by 
hand' another instance of multistationarity in which two nondegenerate steady states exist. 

For the remaining 107 networks, the CRN Toolbox concluded that they do not admit mul- 
tiple steady states. A portion of a report produced by the Toolbox for such a network follows: 

Taken with mass action kinetics, the network CANNOT admit multiple positive 
steady states or a degenerate positive steady state NO MATTER WHAT (POSITIVE) 
VALUES THE RATE CONSTANTS MIGHT HAVE. 

The theoretical underpinning of the Toolbox consists of the Deficiency, Advanced Deficiency, 
and Higher Deficiency Theories developed by Ellison, Feinberg, Horn, Jackson, and Ji (9j [TT] 

Hang. 
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7.3 Classifying irreversible two-reaction networks 



Next, we consider the irreversible versions of the reversible two- reaction networks studied. That 
is, we are interested in the networks obtained from the 386 reversible networks by making one 
or both of the non-flow reactions irreversible (each reversible reaction can be made irreversible 
in two ways). So, each reversible network has 8 relevant subnetworks. Recall that smaller 
sub-CFSTRs, those containing only one directed non-flow reaction or one pair of reversible 



non-flow reactions, were already analyzed in Theorem 2.11, and the bimolecular hypothesis 
ensures that none are multistationary in the setting here. By Theorem 3.1 only subnetworks of 
one of the 35 multistationary reversible networks can be multistationary. Therefore, we must 
examine only 35 * 8 = 280 such networks. We again applied the Toolbox [10] . We found that 
each of the 35 reversible CFSTRs has a unique minimal sub-CFSTR N{ that admits multiple 
positive steady states. Of these 35 subnetworks N{, 29 of them have two directed non-flow 
reactions, while the remaining 6 have non-flow reactions that consist of 1 reversible reaction 
and 1 directed reaction. Examples of both types appear in Figure [2] Thus, a bimolecular 
two-reaction (possibly irreversible) CFSTR admits multiple positive steady states if and only if 
one of these 35 minimal networks is a subnetwork (part 1 of Corollary |7.2[). 



A«A+B 
A+B->2A 



A^A+B 
A+B-> 2 A 




AhA+B 
A+Bm2A 



A^A+B 




A^2A 


A+B«2A 




A<-»2B 



Am 2 A 
Aw2B 



Figure 2: Here we display two of the 35 multistationary bimolecular reversible two- reaction CFSTRs 
(all inflow and outflow reactions are implied), together with all their respective irreversible multista- 
tionary sub-CFSTRs (subnetworks). The poset relation depicted is the subnetwork relation. At left, 
the reversible CFSTR defined by reactions A t» A + B t» 2A has three multistationary sub-CFSTRs, 
which are displayed above. Similarly, the example on the right, defined by 2A ^ A 2B, has only one 
multistationary sub-CFSTR. More generally, each of the 35 such reversible networks admits a unique 
minimal multistationary sub-CFSTR. These minimal subnetworks fall into two classes: 29 of them have 
the form of the example displayed on the left (the minimal network has two directed reactions) , and the 
remaining 6 have the form of the example on the right (the minimal network has one reversible reaction 
and one directed reaction). These 35 minimal sub-CFSTRs appear in Figure [3j 

Finally, we examined the poset obtained from the iVj with respect to the relation of 'embed- 
ded networks' which is displayed in Figure [3} This poset has 11 minimal elements, which are 
the bimolecular two-reaction CFSTR atoms of multistationarity. It follows that a bimolecular 
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Figure 3: Here we display the 35 multistationary bimolecular two-reaction CFSTRs that are minimal 
with respect to the subnetwork relation. The poset relation depicted here is that of embedded networks: 
an arrow points from a network iVtoa network G if N is an embedded network of G. In addition, each 
such edge is labeled by the species that is removed to obtain N from G; for example, C(2) denotes that 
G contains two molecules of species C, and these two are removed from G to obtain N. Two networks 
in the poset are displayed with the same height if they contain the same number of molecules. The 
11 CFSTR atoms of multistationarity are marked in bold/red; they are the networks that have only 
outgoing edges in the figure (at the tops of each component of the poset). All three figures in this work 
were created in Mathematica. 
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two-reaction (possibly irreversible) CFSTR admits multiple positive steady states if and only if 
it contains one of these 11 atoms as an embedded network (part 2 of Corollary 7.2). 

We end by noting that prohibitively many bimolecular t/iree-reaction networks exist, so 
currently there is no classification of those CFSTR atoms. 
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